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Finite size effects in a model for plasticity of amorphous composites 
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We discuss the plastic behavior of an amorphous matrix reinforced by hard particles. A mesoscopic 
depinning-like model accounting for Eshelby elastic interactions is implemented. Only the effect of a 
plastic disorder is considered. Numerical results show a complex size-dependence of the effective flow 
stress of the amorphous composite. In particular the departure from the mixing law shows opposite 
trends associated to the competing effects of the matrix and the reinforcing particles respectively. 
The reinforcing mechanisms and their effects on localization are discussed. Plastic strain is shown 
to gradually concentrate on the weakest band of the system. This correlation of the plastic behavior 
with the material structure is used to design a simple analytical model. The latter nicely captures 
reinforcement size effects in — (log observed numerically. Predictions of the effective flow 

stress accounting for further logarithmic corrections show a very good agreement with numerical 
results. 
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I. INTRODUCTION 

The introduction of foreign particles in an 
amorphous matrix has for long been used for 
the strength reinforcement of disordered mate¬ 
rials [1]. The most classical strategy consists 
in adding rigid particles or fibers in order to 
enhance the elastic properties of the composite 
material. 

An additional or alternative strategy consists 
in the modification of the plastic properties. 
Here, the effect on the overall strength is more 
delicate. The introduction of hard particles in 
a very ductile matrix tends to increase the ef¬ 
fective yield stress, hence the strength. A good 
illustration of this approach can be found in 
the development of materials for road pave¬ 
ments ms]; mineral micrometer scale fillers 
are introduced in a viscous bitumen to make 
it viscoplastic; the obtained mastic asphalt is 
then reinforced through the addition of milli- 
metric to centimetric aggregates. More recently 
the introduction of a ductile phase has been 
used to reinforce metallic glasses 012. In this 
case, the ductility of the second phase enables 
one to control the development of shear-bands, 
thus preventing the nucleation of cracks. A re¬ 
inforcement effect is obtained despite the fact 
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that the effective yield stress of the amorphous 
composite is lowered with respect to that of the 
matrix. 

The understanding of the plastic behavior of 
amorphous composites thus appears to be cru¬ 
cial in the design of modern materials. Efforts 
in theoretical and numerical modelling have 
been recently performed to study the effects 
of micro-alloying in metallic glasses mi] and 
of the addition of aggregates in asphalt mix¬ 
tures [iQHia. 

From the theoretical mechanics point of view, 
the determination of effective mechanical prop¬ 
erties is a matter of homogenization. While 
this field has been intensively explored in case 
of elastic properties [T], results are much more 
scarce for non-linear behaviors like fracture [13 
[T4] or plasticity [T5lll9] . In particular, stan¬ 
dard homogenization approaches fail to account 
for size effects nmiE]. Only the development 
of strain-gradient theories (relying on the in¬ 
troduction of an ad-hoc internal length scale) 
has so far succeeded to reproduce size depen¬ 
dence m- Still, these approaches only pre¬ 
dict the mean behavior and cannot cope with 
sample-to-sample fluctuations. 

Here we develop an alternative approach, 
based on the recent development of depinning- 
like mesoscopic models of amorphous plastic¬ 
ity [2QH25] . The modelling of amorphous plas¬ 
ticity and rheology of complex fluids has seen 
much progress in recent years [23 and a fam¬ 
ily of mesoscopic models EUIMI has emerged 
that rely on two main ingredients: local plastic 
thresholds (amorphous plasticity results from 
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series of local rearrangements of the amorphous 
structure EH EH) and account of elastic in¬ 
teractions (local plastic events occur in a sur¬ 
rounding elastic matrix and induce internal 
stresses m)- These models show scaling prop¬ 
erties close to the effective yield stress (here 
seen as a critical threshold) and thus exhibit 
statistical size effects. Another nice feature of 
these models is their ability to reproduce lo¬ 
calization and shear-banding behaviors [30l |3T] . 
The effect of crystalline inclusions in an amor¬ 
phous matrix has recently been discussed along 
such lines in Ref. [32] with a particular empha¬ 
sis on the localization behavior. 

Here we specialize the model recently pre¬ 
sented in [111 [33] to the case of amorphous com¬ 
posites by considering a binomial distribution 
of local plastic stress thresholds in order to re¬ 
produce the introduction of a fraction of hard 
particles in an amorphous matrix. The simplis¬ 
tic model presented in the following will obvi¬ 
ously not be able to give a realistic account of 
the whole richness of the mechanical behavior 
of amorphous composites. However results are 
expected to be generic for this class of materi¬ 
als. 

In section [HI we introduce the model, in sec¬ 
tion [ITT] we present the complex size depen¬ 
dence of the yield strength measured on the 
amorphous matrix and amorphous composites 
with a growing fraction of particles. In sec¬ 
tion]^ we discuss the hardening mechanisms 
at play in amorphous composites and the lo¬ 
calization behavior. We emphasize in particu¬ 
lar the interplay between the gradual localiza¬ 
tion of the plastic deformation and the build¬ 
ing of a strongly correlated internal stress field. 
Elaborating on the numerical observations, we 
present in section |V| an analytical model that 
accounts quantitatively for the size effects of 
the effective yield strength of amorphous com¬ 
posites. Mathematical details of the model are 
provided in a separated appendix. Our main 


findings are finally summarized in section VI 


II. MODELLING AMORPHOUS 
PLASTICITY: FROM GLASSES TO 
AMORPHOUS COMPOSITES 

The modelling of amorphous plasticity has 
recently given rise to an increasing interest [26] . 
Unlike crystalline plasticity that results from 
the motion of dislocations of the ordered lat¬ 
tice, amorphous plasticity results from series 
of localized rearrangements of the disordered 
structure EHISHIEI- Such local plastic events 






FIG. 1. Effect of a plastic inclusion in an elastic 
matrix submitted to a pure shear biaxial loading 
E = Tiyy = —Exec- Only the spherical inclusion 
experiences plasticity. The first row represents the 
plastic strain the second row the elastic strain 
and the third row the total strain s — . 

The strain fields are represented by a color scale on 
the reference mesh on the left column and on a de¬ 
formed mesh on the right column. One recognizes 
the traditional quadrupolar symmetry associated 
to the Eshelby inclusion. 


induce internal stresses within the surrounding 
material [35l [36] . The latter can be seen as an 
elastic matrix around a plastic inclusion and 
the stress associated to the rearrangement com¬ 
puted in the spirit of the problem of the eigen 
strain early introduced by Eshelby [29] . 

In Eig. we show the effect of such a plas¬ 
tic inclusion on the surrounding matrix. The 
inclusion has experienced a pure shear = 
= ^ 0 ^- The amplitude of the shear strain 
^ = ^xx ~ £yy is represented by a color scale on 
a reference undeformed grid (left column) and 
on a grid deformed according to the total dis¬ 
placement (right column). The first row shows 
the plastic strain non-zero only within the 
inclusion. The second row shows the elastic 
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strain field 5®^ The latter is negative within 
the inclusion. Outside the inclusion it exhibits 
a quadrupolar symmetry: negative along the 
axes at 0 and 90 degrees, and positive along 
the directions at ±45 degrees. The third row 
shows the total strain e = 

The precise expression of the internal stress 
field depends on the details of the plastic 
strain field and the geometry of the rearranging 
region but in the far field, the dominant term 
obeys the universal form: 

= ( 1 ) 

where and A are the mean plastic strain 
experienced by the inclusion and the area of 
the inclusion, respectively. 

The quadrupolar elastic interaction associ¬ 
ated with localized plastic events is the essen¬ 
tial ingredient of the recent models of amor¬ 
phous plasticity and/or rheology of complex flu¬ 
ids. Indeed its anisotropic nature is responsible 
for non-trivial behaviors mEZ]. In particular 
the presence of multiple soft modes of the elas¬ 
tic interaction, corresponding to the existence 
of extended modes of plastic strain which do 
not induce internal stresses, e.g., shear bands, 
is responsible for a complex localization behav¬ 
ior. 


A. A mesoscopic model of amorphous 
plasticity 

Following the model introduced in Refs. [21] 
[24] . the system is discretized on a two- 
dimensional square lattice with a lattice con¬ 
stant that is larger than the typical size of the 
plastic reorganizations. Each site (i,j) has an 
internal stress cr®j, a local plastic threshold cr®^-, 

and a local plastic strain s^j. A pure shear ex¬ 
ternal loading is considered: E®^^ = — E®^^ = 
E®^^. It is assumed that the reorganizations 
at a microscopic scale obey the same symme¬ 
try as the external loading, i.e., a site (i,j) 
undergoes a plastic deformation in pure shear: 
^xx — ~ ^ local criterion of plastic¬ 

ity is considered, the limit of elasticity is thus 
defined for a site (i, j) by: 

+ ( 2 ) 

Values of a® are drawn from a random distri¬ 
bution. No spatial correlations are considered. 

Whenever the criterion is locally satisfied, 
say on site the site undergoes an incre¬ 

mental plastic strain This value is drawn 


from a uniform distribution in [0,do]. To ac¬ 
count for the structural change experienced by 
the rearranging zone, the local plastic thresh¬ 
old is updated to a new value. As discussed 
above, the local plastic event also induces an 
incremental internal stress on every lattice site 

(.i): 

Safj = G^USsl\ (3) 

where the star denotes the convolution opera¬ 
tion and G®^ is a quadrupolar kernel account¬ 
ing for the elastic reaction of the matrix to a 
unit plastic event. Here we consider bi-periodic 
boundary conditions and G®^ is computed from 
Fourier space [24[ |38l [39] . 

The system is driven with an extremal dy¬ 
namics: only one site is deformed per iteration 
step. An iteration step corresponds to (i) iden¬ 
tify the weakest site (ioGo) for a given configu¬ 
ration, (ii) update the plastic strain and 

the plastic threshold at this particular site 
and (Hi) update the internal stress cr®^ all over 
the system. A new configuration is thus ob¬ 
tained and the next iteration can be performed. 
Extremal dynamics [21] is a way of driving the 
system at a vanishing strain rate, very close in 
spirit to the athermal quasi-static driving used 
in some atomistic simulations [40l [41]. Note 
that the same model can be driven with other 
kinds of dynamics, e.g., constant stress, kinetic 
Monte Carlo. 

A direct outcome of a simulation is the evo¬ 
lution of the external stress E®^^ versus the av¬ 
erage plastic strain where the average (•) 

represents the average over the different sites at 
a particular iteration step. The average plastic 
strain is directly proportional to the num¬ 
ber of iteration steps so that {e^^) can be seen 
as a fictitious time. 

In Fig. we give a sketch of a simple plas¬ 
tic behavior. A typical stress-strain curve ob¬ 
tained upon monotonous loading is shown. A 
(reversible) elastic behavior is first observed up 
to the yield stress value E^. Above this value, 
plasticity sets in (a residual plastic strain is ob¬ 
tained upon unloading). The following curva¬ 
ture of the stress-strain curve is characteristic 
of a hardening behavior: if a unloading/loading 
cycle is performed, a new (larger) value of the 
elastic limit is obtained. A stress plateau is 
eventually reached that defines the ultimate 
flow stress E^. 

In the present framework, the external load¬ 
ing is not monotonous. Rather, the external 
stress E®^^ is a fluctuating quantity which is 
adapted at each iteration step so that the crite- 
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FIG. 2. Sketch of a simple plastic behavior. Plas¬ 
ticity sets in at yield stress a hardening stage 
follows until a stress plateau is reached. The latter 
stress value defines the flow stress E^. The plastic 
strain is dehned as the total strain £ minus the 
elastic strain 


rion of the weakest site is satisfied. The macro¬ 
scopic flow stress E^ of a given configuration 
is thus obtained as the maximum value of the 
external stress over the simulation: 

, (4) 

where t is an iteration step. For an external 
loading E®^^ < E^, plastic deformation will 
eventually stop while any loading E®^^ < E^ 
will allow it to develop indefinitely. 


= 1 and = 0.5. The inclusions can be 
either less or more ductile than the amorphous 
matrix. In the cases of interest presented above, 
their nature is often crystalline. We thus as¬ 
sume low fluctuations of the plastic properties 
of the inclusions and we consider that they are 
characterized by a constant plastic threshold, 
(jC = inclusions get the same yield stress 

and this value does not change after an inclu¬ 
sion has experienced plastic deformation. Here 
we restrict the scope to the case of hard par¬ 
ticles: E^ > In order to reduce the space 
of parameters we also consider that the typical 
plastic strain undergone by the inclusions is the 
same as in the amorphous matrix. 

C. Overview of the simulations 

Simulations were performed with sizes rang¬ 
ing from N = 16 up to N = 256, and a num¬ 
ber M = 40 of independent realizations of the 
disorder. The fraction of inclusions was varied 
between, (p = 2.5 x 10“^ and p = 0.99. Dif¬ 
ferent values of the contrast between inclusions 
and matrix were used : E^ = 4,10,40 and the 
value E^ = 10^ was used to mimic infinitely 
hard particles. Most of the following discussion 
will focus on the case E^ = 10. 


III. A SIZE DEPENDENT EFFECTIVE 
YIELD STRESS 


B. Application to amorphous composites 

The model presented above can be easily 
applied to the case of amorphous composites. 
A major hypothesis (already performed in the 
bare model) consists in assuming the homo¬ 
geneity of the elastic properties. Only the effect 
of a plastic disorder will be considered in the 
following. To represent the amorphous com¬ 
posite we consider a fraction p of inclusions 
randomly distributed in an amorphous matrix. 
Here the size of the inclusions is assumed to be 
given by the mesh size so that no correlation is 
considered in the spatial distribution of inclu¬ 
sions. The fraction of inclusions is defined by 
(/> = Vinc/iV^ where Vine is the number of hard 
inclusions and N the linear size of the lattice. 

A bimodal distribution is used to account for 
the respective plastic thresholds of the matrix 
and the inclusions. For the amorphous matrix, 
the plastic threshold is drawn from a uniform 
distribution [a^ — . Here we choose 


A. Amorphous matrix 

We first discuss size effects in the case of a 
mere amorphous matrix, i.e., in the absence of 
hard particles. The ultimate yield strength or 
flow stress E^ of the material is defined as the 
maximum stress experienced by the material 
for a given simulation. 

In Fig. we show the evolution of the ulti¬ 
mate yield strength with the system size. A 
slight decrease is observed. In the inset, we 
show that the evolution is consistent with a sim¬ 
ple power law dependence: 


where E* is the flow stress in the limit of an 
infinitely large system and A is a constant. 
Such a power-law dependence is consistent with 
the depinning-like nature of the model. In 
this context ED EH [37], the plastic flow stress 
can be viewed as a critical threshold between 
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FIG. 3. (color online) Variation of the ultimate 
yield strength with the system size N for a mere 
amorphous matrix (0 = 0) with a yield stress G 
[0.5; 1.5]. The line corresponds to the power law 
expression E^ = E* + ^. As shown in the inset this 
evolution is consistent with the numerical data. 

a static phase (no plasticity) and a dynamic 
phase (plastic flow). The fluctuations of the de- 
pinning threshold measured over a finite length 
scale here simply reflect the divergence of the 
correlation length in the vicinity of a critical 
threshold ^ « 1/ — /*|~^- The present results 
are consistent with the rough estimate u ^ 1 
obtained in previous works [33l[37]- Fig.|^gives 
another illustration of this critical-like behav¬ 
ior. This figure shows the variation of the stan¬ 
dard deviation (5E^ with the average flow stress 
E^. The variation is reasonably reproduced by 
an affine relationship ^E^ = a(E^ — E*). This 
is consistent with the expected critical behavior 
(E* — E^) (X (5E^ oc . The intercept value 

E* can be seen here as the extrapolated value 
of the effective flow stress at infinite size. 

Note that, independently of the system size, 
the values of the effective flow stress lie signif¬ 
icantly above the simple average of the micro¬ 
scopic thresholds = 1. 


B. Amorphous composites 

Second, we discuss the dependence of the ul¬ 
timate yield strength on the fraction of inclu¬ 
sions and on the size of the system. 

Size dependence - In Fig.[^ we show the size 
dependence observed for amorphous compos¬ 
ites with volume fractions of inclusions ranging 
from 0 = 0 to 0 = 0.16. 

For low fractions of hard inclusions, the be¬ 
havior is similar to that obtained for the amor- 



FIG. 4. (color online) Variation of the stan¬ 
dard deviation (5E^ of the ultimate yield strength 
E^ with E^ for the amorphous matrix with a 
yield stress G [0.5; 1.5], for system sizes N = 
16, 32, 64, 128, 256. The standard deviation is 
obtained for 40 realizations. As expected for a crit¬ 
ical transition, a linear behavior is obtained. An 
extrapolation at zero standard deviation gives an 
estimate of the critical threshold, here the yield 
stress E* at inhnite size. 



FIG. 5. (color online) Variation of the ulti¬ 
mate yield strength E^ with the system size N 
for the amorphous matrix (0 = 0) and amorphous 
composites with various fractions of hard inclu¬ 
sions. The yield stress of the amorphous matrix 
is cr^ G [0.5; 1.5], the yield stress of the inclusions is 
E^ = 10. Depending on 0 the yield strength shows 
either a decreasing or an increasing trend with in¬ 
creasing system size. 

phous matrix. The yield strength decreases 
with increasing system size and converges to¬ 
wards a finite value for large system sizes. 

Surprisingly, the behavior is markedly differ¬ 
ent for large fractions of inclusions: the ulti¬ 
mate yield strength increases with increasing 
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system size. At intermediate values of the frac¬ 
tion of inclusions, the evolution of the yield 
strength even appears to be non-monotonic. 

Mixing law -In Fig. we show the evolution 
of with the fraction (j) of inclusions of yield 
stress = 10 for system sizes ranging from 
A" = 16 to A" = 256. The error bars indicate 
the standard deviation computed on the differ¬ 
ent realizations performed for a given pair of 
parameters (0, A). 

A clear size-effect is observed. The curves 
obtained for different values of A do not super¬ 
impose. The larger the system, the larger the 
reinforcement effect induced by the hard inclu¬ 
sions and the closer the effective yield strength 
to the value obtained from a simple linear mix¬ 
ing law: 

A) = (1 - (/))S^(A) + , (6) 

where is the ultimate yield strength of the 
sole amorphous matrix and the yield stress 
of the hard sites. Note that the value ob¬ 
tained from a linear mixing law, known as the 
Voigt average in the context of homogenization 
is usually expected to be an upper bound [T]. 
While this statement is true for homogeniza¬ 
tion of linear properties such as conductivity or 
elasticity, it does not necessarily hold for non¬ 
linear properties such as fracture or plasticity. 
In the latter case, out-of-equilibrium mecha¬ 
nisms may allow the effective property to reach 
values above the Voigt bound [laiii]. 

Although it often fails to reproduce quantita¬ 
tively the experimental data, the simple linear 
mixing law m remains widely used in mate¬ 
rial science to account for the effects of plastic 
reinforcement laiiiiii]. 

Another feature, here emphasized in the inset 
of Fig.|^ can be pointed out: for a given system 
size A, no reinforcement is observed below a 
threshold value (pdN) of the volume fraction of 
hard inclusions. The larger the system size A, 
the smaller the threshold value (pd^). 

Despite its simplicity (scalar model, perfect 
plasticity), the present model is thus character¬ 
ized by a complex behavior. In particular it ex¬ 
hibits a clear size effect that can usually only be 
reproduced in the framework of more complex 
descriptions of plasticity such as strain-gradient 
based theories m- A key ingredient is here the 
account of the elastic interaction induced by the 
local plastic events. 



FIG. 6. (color online) Variation of the ultimate 
yield strength with the fraction p of hard inclu¬ 
sions for a yield stress A G [0.5; 1.5] of the matrix, 
a yield stress E^ = 10 of the hard inclusions and 
for five different system sizes A = 16, A = 32, 
A = 64, A = 128, and A = 256. The same data 
are shown in the inset in semi-logarithmic scale. 


IV. HARDENING AND 
LOCALIZATION 

We now discuss in more details the plastic be¬ 
havior of the model amorphous composites. In 
the following, we try to unveil the mechanisms 
at play in the hardening regime. We shall dis¬ 
criminate between two different effects, respec¬ 
tively associated to a structural evolution of the 
amorphous matrix and a concentration of the 
stresses on the hard particles. We then show a 
gradual localization of the plastic deformation 
on the weakest band of the material. 


A. Stress-Strain curves 

Figure [^displays the stress-strain curves ob¬ 
tained for four different values of the inclusion 
yield stress E^ = 4, 10, 40, 10^ (the latter 
case being meant to mimic infinitely hard in¬ 
clusions) and for different volume fractions p 
ranging from 0 to 0.25. Note that in order to 
emphasize the hardening regime the variation 
of the stress is represented versus the sole plas¬ 
tic strain. 

Two successive hardening regimes can be dis¬ 
tinguished before the stress plateau correspond¬ 
ing to the flow stress is reached. The first one 
is related to the hardening of the amorphous 
matrix. The second one is directly induced by 
the presence of hard particles. 
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FIG. 7. (color online) Stress-strain curves for a 
yield stress of the matrix G [0.5; 1.5], volume 
fractions of hard sites 0 = {0, ...,0.25} and for a 
yield stress of the hard inclusions of = 4 (a), 
E^ = 10 (b), E^ = 40 (c), and E^ = 10® (d). 



FIG. 8. (color online) Evolution of the distribution 
of local plastic thresholds P{cf^) upon plastic de¬ 
formation. A gradual exhaustion effect is observed 
until a stationary distribution is reached. 


B. Statistical hardening of the amorphous 
matrix 

In this subsection, hardening in the pure 
amorphous matrix is considered. At low plastic 
strain, a gradual hardening of the amorphous 
matrix takes place. This phenomenon which 
has been discussed in Refs. mjim results from 
the progressive exhaustion of the weakest sites 
of the matrix. We show in Fig. |^the gradual 
evolution of the distribution of the local plastic 
thresholds P{<J^) upon deformation in the case 
of the sole matrix. The larger the deformation, 
the narrower the distribution and the closer the 
mean to the upper bound value. 

This structural evolution can be understood 
in the following way. After plastic rearrange¬ 
ments, the sites are given a new plastic thresh¬ 
old drawn from the same random distribution 
as the initial ones. The systematic bias be¬ 
tween the weak thresholds of the failing sites 
and the “normal” thresholds that replace them 
after deformation induces an evolutionary-like 
transient increase, reminiscent of self-organized 
criticality models [45]. It can be seen in Fig. 
that at low fractions of inclusions, this exhaus¬ 
tion mechanism is the only one to hold and 
hard particles do not contribute to the rein¬ 
forcement. Indeed the stress-strain curve at 
low volume fractions of hard sites is identical 
to that of the pure amorphous matrix. 

A complementary view of this statistical 
hardening is given in Fig.[^ Here, instead of the 
local plastic thresholds, we show the evolution 
of the distribution of the effective thresholds 
P(cr®^) where = cr® — Indeed, follow¬ 
ing Eq. the local criterion for a given site 
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FIG. 9. (color online) Evolution of the distribution 
of effective plastic thresholds P(crc^) where = 
upon plastic deformation. The sharp lower 
front is associated to the emergence of a global yield 
stress. The latter gradually increases upon plastic 
deformation (hardening) until it reaches its final 
value (the stress plateau of the stress-strain curve). 

(i,jf) can be rewritten as < crF — cr®j. In 
other terms, the local thresholds are dressed by 
the internal stress. Following the evolution of 
the distribution upon deformation, we recover 
the hardening effect. Interestingly, even in the 
transient stage, one can identify a sharp front 
associated to the lower bound of the distribu¬ 
tion. This directly corresponds to the emer¬ 
gence of a yield stress. The disordered system 
has self-organized and in the transient stage one 
can unambiguously define a yield stress that de¬ 
pends on an internal variable, the cumulated 
plastic strain. This also shows the dependence 
of the macroscopic plastic properties on the 
past mechanical history. 


C. Inclusion hardening 

When hard particles are present in the amor¬ 
phous matrix, an additional hardening stage is 
observed. As seen in in Fig.[^ this second stage 
takes place at higher plastic strains than the 
matrix hardening stage. We observe that the 
higher the fraction of hard particles, the sooner 
the onset of the second hardening regime. Note 
also that below a certain fraction of hard inclu¬ 
sions, no second stage of hardening is observed. 

In comparison to the pure amorphous matrix, 
the initial distribution of plastic thresholds is 
bimodal in a composite. In the initial stage of 
the deformation, due to the high contrast of 
plastic thresholds, only sites of the amorphous 
matrix can experience plasticity. The plastic 



FIG. 10. (color online) Evolution of the distri¬ 
bution of effective plastic thresholds where 

(7®^ = cr^ —upon plastic deformation in an amor¬ 
phous composite with 16% of hard particles of yield 
stress = 10. The building of internal stresses 
gradually smears out the peak associated with hard 
particles. Gonversely, the sharp front correspond¬ 
ing to the global yield stress increases upon plas¬ 
tic deformation longer than in the case of the sole 
amorphous matrix. 


events induce internal stresses. Hard particles 
can sustain a level of internal stress much higher 
than that of the amorphous matrix and act here 
as a kind of internal skeleton bearing most of 
the stress exerted on the structure. 

Again it is of interest to follow the distribu¬ 
tion of effective thresholds. In Fig. we show 
the evolution observed for an amorphous com¬ 
posite with 16% of hard particles of yield stress 
E;^ = 10. We see that upon plastic deforma¬ 
tion, the build-up of internal stress on hard par¬ 
ticles has a clear effect: it tends to smear out 
the peak around E^ = 10. In the mean time the 
lower part of the distribution, in particular the 
sharp front that corresponds to the global yield 
stress keeps on increasing. This second hard¬ 
ening stage is much longer than the statistical 
hardening of the amorphous matrix. Stationar- 
ity is eventually obtained when the second peak 
has entirely disappeared. 


D. Localization: no-slip bands 

In order to reveal the hardening mechanisms 
induced by the hard particles, we give a closer 
look at the spatial organization of the plastic 
strain field. In Fig. we show in the top row 
maps of the relative plastic strain 5?^/ {s^^) ob¬ 
tained after a long simulation for three concen¬ 
trations of particles (0 = 10“^, 10“^, 10“^) and 
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FIG. 11. (color online) (a),(b),(c): maps of the relative plastic strain for a system size N — 64, 

a yield stress of hard sites = 10, and volume fractions of hard inclusions 0 = 10“^, 10“^, 10“^ of hard 
sites,respectively. (d),(e),(f): maps of the associated hnal conhgurations of plastic thresholds . The 
positions of hard sites are visible in dark red. 


in the bottom row the associated maps of plas¬ 
tic thresholds (in the final configuration of a 
long simulation) indicating the position of the 
hard sites. 

The low concentration case (panels (a) and 
(d) of Fig. 11, 0 = 10“^) gives a good illus¬ 


tration of the effect of adding hard sites on 
plastic deformation. We see that the plastic 
strain field is not homogeneous. In this exam¬ 
ple where only 3 hard particles are present, we 
observe, as expected, that the hard particles 
are barely deformed. Interestingly, plastic de¬ 
formation is also small along the bands at ±45° 
that intercept the hard sites. Plasticity is inhib¬ 
ited along a set of “no-slip” bands induced by 
the presence of hard particles. These bands ori¬ 
entated at ±45° obviously reflect the symmetry 
of the quadrupolar elastic interaction discussed 
above. While the low fraction of hard inclu¬ 
sions shown in this example is not sufficient to 
induce any reinforcement, it gives a simple clue 
on the strengthening mechanism: hard parti¬ 
cles inhibit the natural slip systems associated 
to the elastic kernel [3^. 


J 


In the medium concentration case (panels 
(b) and (e), 0 = 10“^), the (relative) plastic 
strain field is more heterogeneous. One recov¬ 
ers patterns orientated at ±45° and it is possi¬ 
ble to distinguish between two kinds of bands: 
bands containing hard sites are much less de¬ 
formed than those not containing hard sites. 
In other words, the lattice of no-slip bands is 
much denser and only the sites not intercepted 
by these bands can easily undergo plastic de¬ 
formation. 

In the high concentration case (panels (c) and 
(f), (j) = 10“^), the (relative) plastic strain field 
is highly heterogeneous and actually highly lo¬ 
calized. Most of the plastic deformation con¬ 
centrates onto one single band. This evolution 
is more clearly shown in Fig. where we rep¬ 
resented maps of the incremental plastic strain 
{t ± St) — (t) where St represents 

a few iteration steps such that (5P^)(^t) = 2 and 
(5P^)(t) = 10, 20, 30, 40, 50 and 60. Upon de¬ 
formation, plastic activity appears to become 
more and more localized. 
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FIG. 12. (color online) maps of incremental plastic strain for a volume fraction of hard inclusions 
0 = 10“^, a system size A = 64, a yield stress of hard sites = 10, and different values of the average 
plastic strain 


E. Localization: the weakest band 


We now try to correlate the plastic activity 
with the underlying structure, here represented 
by the landscape of plastic thresholds. As dis¬ 
cussed above, plastic deformation tends to lo¬ 
calize along bands orientated at ±45° that re¬ 
flect the symmetry of the Eshelby quadrupolar 
elastic interaction. Due to statistical fluctua¬ 
tions, not all possible slip systems encounter 
the same number of hard particles. We de¬ 
fine the weakest band and the strongest 

bands S'^max as the bands containing respec¬ 
tively the smallest and the largest amount of 
hard particles among the 2N possible slip sys¬ 
tems. Here we take into account the two possi¬ 
ble orientations. Note again that we consider 
periodic boundary conditions so that all slip 
systems are a priori equivalent. We can now 
compute the fraction of plastic activity occur¬ 
ring in the various bands. In order to high¬ 
light the gradual development of the localiza¬ 


tion, we proceed as in Sec. IV D we consider 
the evolution of the incremental plastic strain 
field ± ^t) — with 5t a few 

iteration steps such that = 2. 


In Fig. we show the evolution with the 
cumulated plastic strain {e^^) of the frac¬ 
tions /min = A5P^(5'Hmin)/(A5P^) aud /max = 
A5P^(5'Hmax)/(A5P^) of the incremental plastic 
strain borne by the weakest and the strongest 
bands respectively for different concentrations 
of hard particles. If the plastic strain field was 
uniformly spread on all bands, one would ex¬ 
pect /min = /max = 1/A (aud UOt 1/2A bc- 
cause of the redundancy between the two pos¬ 
sible orientations at ±45° and —45°). 

In the case of the sole amorphous matrix, 
the weakest band deforms about twice as more 
than the strongest band: /min ~ 2/max- Foi* 
a fraction (p = 10“^, the effect is a bit more 
pronounced but not spectacular yet, we ob¬ 
serve /min ~ 4/max- This ratio remains reason¬ 
ably constant upon deformation. This is con¬ 
sistent with the typical heterogeneity observed 
in Fig. Note that for such a concentration 
the number of particles falls strictly below the 
number of slip systems so that deformation can 
always find a band free of particles to develop. 
No significant reinforcement is expected in this 
case. 

Above some threshold, all slip systems are 
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FIG. 13. Fractions /min and /max of incremental 
plastic strain borne by the weakest and the 

strongest slip systems SB min and SBma^ containing 
the smallest and the largest amount of particles, 
respectively. 


virtually blocked by hard particles. This is the 
case for the two concentrations (p = 0.01 and 
(j) = 0.16 shown in Fig. Here we see a 
dramatic effect: upon deformation, the weakest 
band bears a higher and higher fraction of the 
plastic activity. Eventually most of the plas¬ 
tic strain occurs within this weakest band. We 
thus observe a strong correlation between struc¬ 
ture and plastic behavior: plastic deformation 
gradually concentrate onto the weakest slip sys¬ 
tem, characterized by the smallest amount of 
hard particles. 


V. A SIMPLE ANALYTICAL MODEL 

The mechanism of reinforcement can thus be 
understood in the following way. Hard particles 
inhibit slip systems. No reinforcement occur 
until all slip systems are blocked. Above the 
associated threshold concentration, all slip sys¬ 
tems are hindered by hard particles and plas¬ 
tic strain gradually localizes onto the weakest 
one, i.e., the one that contains the fewest hard 
particles. The macroscopic plastic behavior is 
thus controlled by the properties of this weakest 
band. In the following, we discuss these two as¬ 
pects, elaborate a simple analytical model and 
compare its prediction with our simulations. 
Mathematical details are presented in a sepa¬ 
rated appendix. 


A. Percolation 

As discussed above, no reinforcement is ex¬ 
pected until all slip systems are blocked by at 
least one particle. Here the two families of 
slip systems associated to the two directions 
at ±45° should a priori be considered. For 
the sake of simplicity, we consider in the fol¬ 
lowing only one of the two orientations. This 
approximation allows to recover a simple one¬ 
dimensional percolation problem. 

We assume here that the distribution of par¬ 
ticles is not spatially correlated and take the 
volume fraction p as the probability for one in¬ 
clusion to be hard. The probability to have ex¬ 
actly n hard inclusions in one randomly chosen 
diagonal is then 

P{Na =n)= r(l - , (7) 

where is the random variable counting the 
number of hard sites on a diagonal, = n is 
the event “n hard sites on a diagonal”, N is 
the number of sites in a diagonal, which is ex¬ 
actly the system size in the square lattice con¬ 
sidered here. We recognize a binomial distribu¬ 
tion. The probability of having at least 1 hard 
inclusion on a diagonal is: 

P(Na >l) = l-(l-c/>)^ . (8) 

There are N diagonals with the same orienta¬ 
tion. They are independent. Consequently, the 
probability to have at least 1 hard inclusion on 
each diagonal is 

p(i?) = (9) 

where the letter B stands for “blocked”. This 
probability is the equivalent of the probability 
of percolation. It is plotted for different system 
sizes versus the volume fraction of hard inclu¬ 
sions in Fig. The probability of having at 
least one hard inclusion per diagonal increases 
with the volume fraction of hard inclusions un¬ 
til it reaches 1. 

The threshold for percolation, or here for all 
diagonals to be blocked, is the volume frac¬ 
tion (j) for which the probability P{B) is the 
steepest. In other words, the threshold for the 
transition corresponds to the volume fraction 
of hard inclusions for which the second deriva¬ 
tive of P{B) vanishes. This volume fraction is 
called critical and denoted 0c- It is equal to 

4>c{N) = 1 - 
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FIG. 14. (color online) Variation of the probability 
P{B) of having at least one hard inclusion on each 
diagonal, dehned in Eq.[^ with the volume fraction 
0 of hard inclusions for different system sizes N — 
16, V = 32, V = 64, V = 128, N = 1024. Inset: 
Variation of the critical fraction (j)c defined in Eq. |10| 
with the system size N. 

The inset of Fig. shows the variation of 
the critical fraction (pc with the system size 
N. The thresholding effect is nicely illustrated 
in Figure where we show the rescaled ul¬ 
timate yield strength (TF{(p,N) — — 

/(pc{N) versus the rescaled volume fraction 
(j)/(j)c{N)^ for different system sizes. This plot is 
to be compared with the inset of Fig. [^show¬ 
ing the same quantities without the rescaling 
by (j)c{N). In the main plot, the curves corre¬ 
sponding to different system sizes collapse onto 
a single master curve, showing that our inter¬ 
pretation of the transition is valid. 


B. Yield stress of the weakest band 

1. Restriction to elastic line depinning 

In section |rV| plastic deformation was shown 
to concentrate onto one single band, the one 
containing the smallest amount of hard parti¬ 
cles. It is thus natural to use the ultimate yield 
strength of that weakest band as an estimate 
for the ultimate yield strength of the whole 
amorphous composite. Ignoring the residual 
plastic strain undergone outside the band, the 
problem is thus reduced to a one-dimensional 
elastic depinning problem very similar to the 
propagation of a crack front in a random land¬ 
scape [IS146]. 

Indeed, if we denote {zi) the plas¬ 

tic strain in the weakest band at location Zi 
where 2 : is the distance along the band, any lo¬ 



FIG. 15. (color online) Variation of the rescaled 
yield strength (E^(0, V) - E^)/(E^ - E^)/0c(V) 
with the rescaled volume fraction p/pc^N) for a 
yield stress E^ = 10 of hard inclusions and differ¬ 
ent system sizes N = 16, 32, 64, 128, 256. Inset: 
Variation of (E^(0, N) - E^)/(E^ - E^) with p for 
the same systems. 

cal plastic strain increment induces along 
the band an internal stress proportional to an 
elastic kernel which is nothing but the restric¬ 
tion on a diagonal of the Eshelby quadrupolar 
stress defined in Eq. More specifically the in¬ 
ternal stress at location Zj induced by the plas¬ 
tic strain increment at location Zi amounts to: 

8cj^, = -A^e^ , 

where Aq and Ai are positive constants. One 
recognizes here the elastic interaction associ¬ 
ated to the trapping of an interfacial crack 
front [T7] . 

The determination of the effective toughness 
of an interfacial crack propagating in a random 
landscape (which also amounts to the critical 
threshold of a long-range elastic line) has re¬ 
cently been discussed in Ref [14]. While the 
effective toughness can significantly exceed the 
simple arithmetic average of the microscopic 
properties when the disorder is highly fluctuat¬ 
ing in the direction of propagation (strong pin¬ 
ning), a simple mixing law is recovered when 
the microscopic toughness is only slowly vary¬ 
ing along the direction of propagation (weak 
pinning). 

In the present case, the hard sites are persis¬ 
tent, i.e., the value of their yield stress does not 
change upon deformation. Besides, the fluctu¬ 
ations of the local thresholds that characterize 
the amorphous matrix are weak compared with 
the yield stress of the hard sites. Weak pinning 


























13 


conditions can thus be considered and a simple 
mixing law used to compute the effective yield 
stress of the band. 


2. How weak is the weakest hand? 


The effective yield stress of the weakest band 
CTy" is thus simply written: 




1 — m_A m_TT 

-H- 

N N 


( 12 ) 


where is the yield stress of the amorphous 
matrix and that of the hard sites and m 
is the number of hard sites in the band. The 
estimate of the ultimate yield strength of 
the material is given by the ensemble average: 




1 - {m) I {m) 
N N 


(13) 


where (m) is the average minimum number of 
hard sites on a diagonal of size N for a fraction 
(j) of hard sites. In the following we define / = 
(m) /N^ the effective fraction of hard sites in the 
weakest band. As it immediately appears from 
Eq.[^ within the weakest band approximation, 
the difference between the effective fiow stress 
and the mixing law estimate directly 
stems from the difference between / and 0. 

The distribution of the number m of hard 
sites is given by the binomial distribution of pa¬ 
rameters (j) and N. An exact formula for the av¬ 
erage (m) of the minimal number of hard sites 
on a diagonal when N diagonals are considered 
is given in the appendix. However, this formula 
contains an infinite sum and is not easy to han¬ 
dle. In order to estimate this minimal value we 
shall resort to an argument of extremal statis¬ 
tics. Depending on the value of 0, the binomial 
converges at large N either to a Gaussian or to 
a Poisson distribution. In the present case we 
are interested in the large deviations of the bi¬ 
nomial distribution [48]. We use recent results 
on the general approximation of the binomial 
distribution gaisni obtained in the context of 
cryptology studies: 


P{m < fN) = 


f r,-ND{f\\4>) 

{^-f)VWf 


(14) 

for A" ^ oc where D(/||0) is the Kullback- 
Leibler divergence defined as: 


^(/Il0) =/lnI + (l-/)lnI—I . (15) 


Here, the fraction / of hard sites in the weakest 
band is estimated via a simple extremal statis¬ 
tics argument: 

P{m < /IV) « 1 . (16) 

Detailed calculations based on the asymptotic 
expansions given in Ref. [50| are presented in 
the appendix. They allow us to obtain an es¬ 
timate of the distance between the fraction / 
(the fraction of hard sites in the weakest band) 
and the parameter (j) of the binomial distribu¬ 
tion (the mean fraction of hard sites): 


f = 4> 


20(1-(/))^ N , , , 

log^(l + ^iv), (17) 


N 




where 

^ I log(2hAr) 
2 2hj\[ “h I 

hN = log . 

V 27r 


(18) 


This immediately sets the distance of the fiow 
stress to the mixing law value E^(0, N) ob¬ 
tained by Eq. 


E^c/., N) = E^(/, - f) [eH - EA(iV)] . 

(19) 

In particular we obtain a clear size effect: 
the convergence to the mixing law scales as 
(logA/A)^/^. This result is illustrated in 
Eig. where we display the variation of the 
rescaled fiow stress cr/^(0, N) with (log A/A)^/^ 
for various values of the fraction 0 of hard sites. 
The rescaled fiow stress (7/^(0, A) is defined as 
the reinforcement factor with respect to the 
fiow stress of the amorphous matrix: 




E^(0, A) -EA(A) 
EH-EA(A) 


( 20 ) 


In the framework of our approximation, we ex¬ 
pect (7/^(0, A) = /((0, A). In particular, fol¬ 
lowing Eq. [l7| we should recover 0—(7/^(0, A) ex 
(logA/A)^7A As shown in Eig. 16, this scal¬ 


ing is nicely obeyed for moderate values of 0. 
Departures from the predicted scaling behavior 
become significant at low values of 0 and A, 
because the analytical estimation holds only in 
the limit of large A and intermediate values of 
0. A numerical estimation of the average num¬ 
ber (m) of hard sites in the weakest band is 
discussed in the appendix and shows that the 
approximation holds surprisingly well even for 
low values of 0 and small system sizes. 
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FIG. 16. (color online) Size scaling of the rescaled 
flow stress aR{(j),N), defined in Eq. of amor¬ 
phous composites with concentration of hard par¬ 
ticles ranging from cj) = 0.04 to cj) = 0.5. Filled 
symbols on the vertical axis correspond to the in¬ 
finite size limit, i.e., the result of the mixing law. 
Indicative dashed lines show the expected scaling 
behavior in {logN/NY^"^. 



FIG. 17. (color online) Effect of the concentra¬ 
tion of hard particles on the rescaled flow stress 
of amorphous composites for different system sizes 
N = 16, 32, 64, 128, 256. The straight line corre¬ 
sponds to the mixing law expected at inhnite size. 
The dashed lines are the analytical predictions of 
Eq. accounting for logarithmic corrections. 


Beyond the prediction of the scaling behav¬ 
ior, the logarithmic corrections accounted for 
in Eq. allow us to test quantitatively our 
predictions for the reinforcement effect of hard 
inclusions in an amorphous matrix. In Fig. ^7\ 
we compare analytical predictions and simula¬ 
tion results for the rescaled flow stress (7/^(0, N) 
with respect to the fraction of hard sites 0. 
Again, our prediction of effective flow stress 
happens to be very precise for moderate values 
of (j) and large system sizes. 


VI. CONCLUSION 

The plastic behavior as described in the 
mesoscopic simulations shows two types of sys¬ 
tem size dependence. The first type corre¬ 
sponds to an effect of the amorphous matrix 
only and results from the critical character 
of the yielding transition. In this case, the 
ultimate yield strength decreases with an in¬ 
creasing system size, as 1/N. This system 
size dependence has already been addressed in 
Refs. [33l[37]- A similar critical behavior has re¬ 
cently been advocated in the related framework 
of compressive strength of brittle heterogeneous 
materials m 

The second type of size effect is specific to 
the composite material. Below a critical vol¬ 
ume fraction of hard inclusions depending on 
the system size, no hardening behavior of the 
second type is observed. Above this critical vol¬ 
ume fraction, the hardening behavior depends 
on the system size: the ultimate yield strength 
increases with an increasing system size, as 
— (log We showed that the threshold¬ 

ing effect observed in the simulations is close 
to a percolation transition. We also showed 
that during this second hardening regime, most 
of the plastic strain is concentrated onto the 
weakest band. Therefore, we proposed a sim¬ 
ple model to describe the dependence of the 
ultimate yield strength on the system size 
and the volume fraction. This model is based 
on the assumption that the weakest band bears 
all the plastic strain and governs the ultimate 
yield strength of the entire system. The 
ultimate yield strength is then given by a 
combination of the yield strength of the pure 
matrix and of the hard inclusions, weighted re¬ 
spectively by the fraction of matrix sites and 
of hard inclusions in the weakest band. Us¬ 
ing extremal statistics arguments, we proposed 
an analytical estimate of the average number of 
hard inclusions in the weakest band in the limit 
of large system sizes. The comparison of the an¬ 
alytical estimate with the simulation results is 
satisfactory and our model consequently makes 
a direct link between the structure, represented 
by the plastic threshold on each site, and the 
mechanical behavior. 
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of the average minimum of hard sites on a di- agonal. 
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APPENDIX: ANALYTICAL APPROACH 

1. Exact formula 


To obtain an exact formula for the average 
minimum number (m) of inclusions per diago¬ 
nal, we use its definition: 


N 


(m) = nP{m = n) , 

(21) 

n=0 


N 


= ^ nP{m = n) . 

n=l 

(22) 


The variable n in this expression can be refor¬ 
mulated as a sum: 

N n—1 

n=l k=0 
N-1 N 


= E E Pi^ = n), (24) 

k=0 n=k-\-l 

(25) 


where the indices in the two sums are enumer¬ 
ated in two different and equivalent ways. The 
sum over the index n can then be contracted 
into: 


N-l 


P{m > /c + 1) , 

k=0 

N 

(26) 

^P{m > n) , 

- 1 

(27) 

^ N 

E(w>«)) , 

n=l 

(28) 

^ N 

^(/^(n,7V-n+l)) , 

(29) 


n=l 


2. Analytical estimation 

To obtain an analytical estimate of the av¬ 
erage minimum number (m) of inclusions per 
diagonal, we use a result from extreme value 
theory on the minimum of N independent and 
identically distributed random variables: 

P{Nd<{m)) = j^, (30) 

where is the random variable counting the 
number of inclusions in any diagonal. In our 
case, the diagonal have N sites. We introduce 
the ratio / = {m)/N convenient in the limit of 
large N and Eq. [3Q| becomes: 

P{Na<fN) = j^, iV^cx). (31) 

As the random variable has a binomial 
distribution with parameters N and 0, we 
can employ a result of cryptography given in 
Refs. j4^l5Q] : 


P{Nd < fN) 


~ _ ^-ND(f\\4>) 

{(p- f)^/WNf 


(32) 

for A ^ oc where D(/||0) is the Kullback- 
Leibler divergence defined as: 


il(/||0)=/lnl + (l-/)ln^. (33) 


We then introduce 6 = 0 — /; and e is expected 
to tend to 0 in the limit of large A, i.e., / tends 
to (/) in the limit of large A. Using Eqs. 
and ^and another result from Refs. @21EO] on 
the behavior of the Kullback-Leibler divergence 
valid for 0(e) = 0{e/(j)) = 0(e/(l — (j))): 

^(0-e||<(>) = 2^(f_^) +Q(^')^ (34) 

we obtain: 


— = ^ - e 2c/>(i- 

N ei/27rA(0 — e) 

To first order in e, this becomes: 


(35) 


where is the random variable counting the 
number of hard sites in a diagonal and 1(f) the 
regularized incomplete beta function. The reg¬ 
ularized incomplete beta function is used here 
to express the cumulative distribution of the 
binomial distribution. The formula Eq. ^9 


IS 


explicit but hard to evaluate for large values of 

A. 


1 = 

A e 


(36) 


We define e' = e/y^0(l — 0)/A and Eq. 
equivalent to 


36 


IS 


= 21n- 


A 


\/^e' 


(37) 
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To obtain an approximate solution to this tran¬ 
scendental equation, we define the variable r 
such that: 


e' = \/2/iiv(l + r) , 


(38) 


where = \n -^=. The variable r tends to 0 

V 27r 

in the limit of large N. We also have: 

e'2 = 2hN{l+rf = 2hN{l + 2r) + 0{r^) , (39) 

to first order in r. Using the transcendental 
equationand iterating once in e, we get: 

N 


= 21n 


\/^\/2hiv(l + r) 


(40) 


Equating the right hand sides of Eqs.[39|and 40 
leads: 


2hAr(l + 2r) = 2 In 


= In 


N 


(41) 


\/^\/2hAr(l + r) 

N 

21n(l + r) , 


= In 


■\/27r-\/2hN 

N 

-s/27r \/2h jsf 


- 2r . 


(42) 

(43) 


to first order in r. Eq. is a linear equation 
in r, its solution reads: 


In 


N 


r = 


\/^y/2h]Sf 


— h 


'N 


2h]\[ T 1 


(44) 


According to Eq. it is equal to 


g{N) 


0(1 - 0 ) 

N 


(47) 


Eig. \T8\ displays the function g as obtained for 
10,000 numerical iterations of N drawings from 
a binomial distribution with parameters N and 
0, for different values of N and 0. The analyt¬ 
ical estimation of g{N) given in Eq. 47 works 
very well down to A" = 16 for values of 0 in 
0.1 < d) < 0.9. It is expected as we used Eq.|34 



FIG. 18. (color online) Variation of the function 
g{N) = (0-/)/(l + r(A))/V^with Afor differ¬ 
ent values of 0. The opened squares correspond to 
a numerical evaluation of g{N) using 10, 000 iter¬ 
ations of A drawings from a binomial distribution 
with parameters A and cj). The solid and opened 
squares of the same color correspond to results for 
(j) and (1 — 0), respectively. The solid lines corre¬ 
spond to the analytical estimation given in Eq. 
in the limit of large A. 


Einally, the ratio / = {m)/N can be expressed 
in terms of r(A): 


/ = <(>- y+ r{N)) . (45) 

The solution given in Eq. can be checked 
numerically. The function g is defined by: 


g{N) 


(!>-f 1 

1 H- r{N) y/2h^ 


(46) 


valid for 0(e) = 0{el(j)) = 0{e /(1 — 0)). Eor ex¬ 
treme values of 0 the analytical estimation only 
gives satisfactory results for very large A. The 
0-(l — 0) symmetry is then recovered. How¬ 
ever, the numerical evaluation of the minimum 
is in very good agreement with the simulation 
results for all values of 0 as is shown in Eig.p!^ 
In this figure, the fraction / of hard sites in the 
weakest band is estimated from the simulation 
results as (E^(0, A) - EA(A))/(E^ - ^^(A)). 
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FIG. 19. (color online) Variation of the function 
g{N) = (0 — /)/(! + r{N))/\/2hN with N for dif¬ 
ferent values of 0. The solid circles correspond to 
simulation results. The opened squares correspond 
to a numerical evaluation of g{N) using 10, 000 iter¬ 
ations of N drawings from a binomial distribution 
with parameters N and cj). The solid lines corre¬ 
spond to the analytical estimation given in Eq. 
in the limit of large N. 







